//==============================================================================
// SequenceAlignment.hh
//
//==============================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8
//

#ifndef _SEQUENCE_ALIGNMENT
#define _SEQUENCE_ALIGNMENT

#include "LocalSTL.hh"
#include "Table.hh"
#include "Sequence.hh"

#define MATCH_SCORE  1
#define MISMATCH_SCORE  -1
#define GAP_SCORE  -2
#define GAP '-'

const int TOP=0x00000001;
const int LEFT=0x00000010;
const int DIAG=0x00000100;

//-----------------------------------------------------------------------------
// * SequenceAlignment
//-----------------------------------------------------------------------------
//
class SequenceAlignment{

public:

	SequenceAlignment();
	SequenceAlignment(Sequence S,Sequence T);
  	SequenceAlignment(string S, string T);
  ~SequenceAlignment();


  Sequence mS,mT;
  Sequence mAlignS,mAlignT;
  
  int mRows,mCols;
  int mBRows, mBCols;
  int mIgnorN;

  // !!! KJD For speed, these should be implemented as a valarray!!!
  Table<int>  mScoreTable;
  Table<int>  mPtrTable;
  Table<int>  mbandScoreTable;
  Table<int>  mbandPtrTable;

  int 		  mMaxScore;  //store the best score found
  int 		  mMaxRow;    //alignemnt end point
  int 		  mMaxCol;    //alignment end point
  
  int mAlignStart; // Local alignment starting position for query
  int mAlignEnd;  //Alignment end position for query
  int mTAlignStart; //Alignment end pos for target
  int mTAlignEnd;   //Alignment end pos for target
  int mTBandStart;
  int mTBandEnd;
  

  //match stats, the length of match can be accessed by calling mAlignS.size()
  int mMismatchCount;  //number of mismatch found
  int mSGapCount;       //number of gap in Sbjct seq
  int mTGapCount;      //number of gap in target seq

  // Score function parameters. Default set 
  int mMatchScore;
  int mMismatchScore;
  int mGapScore;
  int mGapExtensionScore;

  // band width
  int mbandwidth;
  int mBandScore;  //the minimal score that the band will extend

  // Score Function
  void	SetScoreParams(int match,int mismatch, int gap);
  void	SetScoreParams(int match,int mismatch, int gap,int extension);
  
  // set if treat N as any base
  void SetIgnorN(int choice);
  int IgnorN();
  //set the minimal band score
  void SetBandScore(int score);

  // Global Alignment
  void 	GetGlobalAlignment(Sequence &AlignS,Sequence &AlignT,int &score);
  void 	ComputeGlobalAlignmentScore();


  // Local Alignment
  int   GetLocalAlignment(Sequence &Source,Sequence &Target,int &score);
  void	ComputeLocalAlignmentScore();  // Whichever is in use now

  // Internal Alignment
  void GetInternalAlignment(Sequence &Source,Sequence &Target);
  void ComputeInternalAlignmentScore();
  void RecoverInternalAlignment();
  void GetBandInternalAlignment(Sequence &Source, Sequence &Target, int targetstart, int bandsize);
  int ComputeBandInternalAlignmentScore();
  void RecoverBandInternalAlignment();

  //end Alignment
  int GetEndAlignment(Sequence &Source,Sequence &Target);
  void ComputeEndAlignmentScore();
  int  RecoverEndAlignment();
  int GetBandEndAlignment(Sequence &Source, Sequence &Target, int targetstart, int bandsize);
  int  ComputeBandEndAlignmentScore();
  int RecoverBandEndAlignment();
  
  //pure overlap alignment
  int GetBandOverlapAlignment(Sequence &Source, Sequence &Target, int targetstart, int bandsize);

  void	PrintPtrTable();
  void	PrintScoreTable();
  void  PrintBandScoreTable();
  void	Print(int start,int end);
  void	Print();
  void  PrintAlign(ostream& output);

  //protected:

  // Generic Utility (KJD Should be moved elsewhere)
  inline int Max(int a,int b,int c);
  inline int Max(int a,int b,int c,int d);
  
  // Local Utility
  void   Swap(int &a,int &b);
  void	PrintSpacePaddedNum(int x);
  int 	Score(char A, char B);
  //inline int FindMaxScore(int i,int j,char*S,char*T);
  //int 	Score(char A,char B);
  int FindMaxScore(int i,int j,char*S,char*T);
  void 	RecoverGlobalAlignment();
  void 	RecoverLocalAlignment();
  void    RecoverLSpaceGAlignment();


  inline void	Set(int signal,int &c);
  inline void	Set(int signal,char &c);
  void	Clear(int signal,int &c);
  inline bool	IsSet(int signal,int c);
  bool	IsSet(int signal,char c);

  void StartTiming(int &Start);
  int  EndTiming(int Start);

};//-----| End SequenceAlignment Class |------------------------------------

//-----------------------------------------------------------------------------
// * LSGAffineGapAlignment
//-----------------------------------------------------------------------------
// Linear Space Global alignment
//
class LSGAlignment : public SequenceAlignment{
public:

  LSGAlignment(){};
  LSGAlignment(Sequence S,Sequence T):SequenceAlignment(S,T){};

  	LSGAlignment(string S, string T):SequenceAlignment(S,T){};
	~LSGAlignment(){};


  // LSpaceGAlignment
  void 	GetAlignment();
  void   HandleLimitingCase(int x1,int y1,int x2,int y2);
  void   GetAlignment(int x1,int y1,int x2,int y2);
  int    FindInterceptWithBestPath(int u,int x1,int y1,int x2,int y2);
};

//-----------------------------------------------------------------------------
// * LSLAffineGapAlignment
//-----------------------------------------------------------------------------
// Linear Space Local alignment
//
class LSLAlignment : public SequenceAlignment{
public:

	LSLAlignment(){};
	LSLAlignment(Sequence S,Sequence T):SequenceAlignment(S,T){};
  	LSLAlignment(string S, string T):SequenceAlignment(S,T){};
	~LSLAlignment(){};

  void 	GetAlignment();
  void   HandleLimitingCase(int x1,int y1,int x2,int y2);
  void   GetAlignment(int x1,int y1,int x2,int y2);
  int    FindInterceptWithBestPath(int u,int x1,int y1,int x2,int y2);
};


// What I really need is LSLocalAffineAlignment

//-----------------------------------------------------------------------------
// * LSGAffineAlignment
//-----------------------------------------------------------------------------
// Linear Space Global Affine gap penalty Alignment. 
//
class LSGAffineAlignment : public SequenceAlignment{
public:
  
  LSGAffineAlignment(){
	 SCurrent = 1;SPrevious = 0;    // score table columns		
	 IxCurrent = 2;IxPrevious = 3;  // Insert X table columns	
	 IyCurrent = 4;IyPrevious = 5;  // Insert Y table columns
	 CsCurrent = 6;CsPrevious = 7;    // "c" table columns
	 CxCurrent = 8;CxPrevious = 9;    // "c" table columns
	 CyCurrent = 10;CyPrevious = 11;    // "c" table columns
  };

  // Pointers for columns
  int SCurrent,SPrevious;    // score table columns		
  int IxCurrent,IxPrevious;  // Insert X table columns	
  int IyCurrent,IyPrevious;  // Insert Y table columns
  int CsCurrent,CsPrevious;    // "c" table columns
  int CxCurrent,CxPrevious;    // "c" table columns
  int CyCurrent,CyPrevious;    // "c" table columns

  // AffineLSpaceGAlignment
  //void 	GetAffineLSpaceGAlignment();
  //void   AffineLSpaceGAlignment(int x1,int y1,int x2,int y2);
  //int    AffineFindInterceptWithBestPath(int u,int x1,int y1,int x2,int y2);

  //void   AffineLSpaceGAlignmentL(char a,int x1,int y1,int x2,int y2);
  //void   AffineLSpaceGAlignmentR(char a,int x1,int y1,int x2,int y2);
  //int    AffineFindInterceptWithBestPath2(char s,char & d,int u,int x1,
  //														int y1,int x2,int y2);

  void GetAlignment();
  void Align(char scaseTL,char scaseLR, int x1,int y1,int x2,int y2);
  void InitScoreTableColumns(int rows);
  void InitColumnTops(int xIdx);
  void FillWithBestScore(int xIdx,int yIdx,int &FromWhere);
  void FillCTables(int yIdx,int CsOrigin,int CxOrigin,int CyOrigin);
  void SwapColumns();
  int GetTraceBack(int rows,int CsOrigin,int CxOrigin,int CyOrigin,int &dir);
  int FindInterceptWithBestPath(char scase,
										  char BestScoreOrigin,
										  int u,int x1,
										  int y1,int x2,
										  int y2);
  

};

#endif


